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We explore bright soliton solutions of ultracold Bose-Fermi gases, showing that the presence of 
p-wave interactions can remove the usual collapse instability and support stable soliton solutions 
that are global energy minima. A variational model that incorporates the relevant s- and p-wave 
interactions in the system is established analytically and solved numerically to probe the dependen- 
cies of the solitons on key experimental parameters. Under attractive s-wave interactions, bright 
solitons exist only as meta-stable states susceptible to collapse. Remarkably, the presence of repul- 
sive p-wave interactions alleviates this collapse instability. This dramatically widens the range of 
experimentally-achievable soliton solutions and indicates greatly enhanced robustness. While we fo- 
cus specifically on the boson-fermion pairing of *^Rb and '*''K, the stabilization inferred by repulsive 
p-wave interactions should apply to the wider remit of ultracold Bose-Fermi mixtures. 

PACS numbers: 03.75.Lm, 67.85.Pq 



I. INTRODUCTION 

The study of wave mechanics and propagation in non- 
linear media is a fundamental concept within physics. In 
particular, solitons are a general non-dispersive solution 
to the one-dimensional non-linear wave equation. Bright 
solitons have been observed in many areas of physics, 
such as in water [T], liquid hydrogen [2], optics [3] and 
atomic Bose-Einstein condensates (BECs) [IHS]. In the 
latter case, these self-trapped matter waves are supported 
by a balance between attractive atomic scattering inter- 
actions and repulsive zero-point kinetic energy between 
bosons. As well as being of fundamental interest in many- 
body quantum physics, bright matter-wave solitons are 
being touted for potential uses in atom interferometry 
[H El H Ej and surface characterization [3]. However, 
when realized in three-dimensions these solitons exist 
only as meta-stable states prone to catastrophic collapse 

Ultracold Bose-Fermi (BF) gases have received a great 
deal of recent experimental attention and have been re- 
alised through ^Li-^Li [12, 23Na-6Li [Tg, ^^Rb-^^K [Hj, 
i74Yb_i73Yb [16] and ^^Sr-^^Sr Tf mixtures. At such 
low temperatures, the scattering of atoms with non-zero 
relative angular momentum is heavily restricted such 
that p-wave and higher interactions are typically negli- 
gible. Furthermore, the Pauli exclusion principle forbids 
identical fermions from interacting via s-wave collisions. 
Thus, for an ultracold Bose-Fermi mixture (in which the 
fermions are identical), the dominant interactions are s- 
wave boson-boson and boson-fermion interactions. It has 
been shown theoretically that a repulsive Bose gas and a 
non-interacting Fermi gas co-existing in a radial waveg- 
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uide can be coupled together by an attractive boson- 
fermion interaction to form a self-trapped state [TU [19] . 
It is these 'Bose-Fermi' solitons that we will consider in 
this paper. We note that the ^^Rb-^°K system appears 
particularly well-suited to support BF solitons since its 
natural boson-fermion interaction is strongly attractive. 
However, just as in the case of an attractively-interacting 
3D Bose-Einstein condensate [TOl dH 120]) a 3D Bose- 
Fermi system is prone to collapse when the interspecies 
interaction becomes too attractive [T51[22HMj . While cer- 
tain effects have been highlighted to raise the threshold 
for collapse, e.g. finite temperature effects in Bose-Fermi 
mixtures |25j and Feshbach resonance management [26] . 
possession of angular momentum [27] and fragmentation 
[55] in Bose-Einstein condensates, the instability to col- 
lapse ultimately remains in the system. 

By exploiting scattering resonances, the p-wave inter- 
action between fermions can now be experimentally en- 
gineered to significant values [H] . It is thus the rationale 
of this work to explore the way in which p-wave inter- 
actions may modify Bose-Fermi soliton solutions. Our 
framework to perform this is a variational approach using 
a gaussian ansatz for the boson and fermion density dis- 
tributions. Via the Gross-Pitaveskii model for the bosons 
and the Thomas-Fermi approximation for the fermions, 
we derive an analytic form for the energy of this cou- 
pled system up to p-wave interactions (for the boson- 
fermion and fermion-fermion interaction) and from this 
we obtain the variational soliton solutions. First we re- 
view the properties of the soliton solutions in the ab- 
sence of p-wave interactions, in agreement with previous 
studies [m [19] , before considering the modifications in 
the presence of p-wave interactions. Specifically we focus 
our results on a ®'^Rb-''''K mixture, due to the naturally 
large and attractive s-wave scattering length between the 
species [TU [151 E] , and fermion-fermion p-wave interac- 
tions, due to their capacity to be engineered experimen- 
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tally [25] ■ 

II. VARIATIONAL MODEL OF BOSE-FERMI 
SOLITONS 

A. System overview and the variational ansatz 

We consider a degenerate gas of identical fermions co- 
existing with a Bose-Einstein condensate of bosons, all 
at zero temperature. Neglecting quantum and thermal 
fluctuations, we will model the fermion and boson gases 
within the mean-field picture. Each gas is confined by 
an axially-homogeneous waveguide potential T^efF} (r) — 
imB{p}Ci;||pjr^, where Wb{f} is the radial trap frequency 
experienced by the bosons {fermions} and TOb{f} is the 
boson {fermion} mass. Due to the low energy of the 
atomic collisions, the s-wave and p-wave interactions are 
modelled by contact interactions characterised by a sin- 
gle length scale, the scattering length. Within the Bose 
gas, the atoms interact predominantly via s-wave scatter- 
ing with characteristic length ob (p-wave interactions are 
negligible). Within the Fermi gas, s-wave interactions are 
suppressed via the Pauli exclusion principle and the lead- 
ing atomic interaction is p-wave with a scattering length 
Of [23]. For overlapping clouds, the bosons and fermions 
additionally interact with each other, predominantly via 
the s-wave interaction, of lengthscale cbf^, but we will 
also include the corresponding p-wave interaction, with 
effective scattering length aBPp |23]. 

Bright solitons require an attractive interaction to en- 
able self-trapping of the wave. Here we shall consider 
the case where this interaction arises from the s-wave 
boson-fermion coupling. A rudimentary requirement is 
thus that the boson and fermion gases are overlapping 
in space and this enables us to assume the same ansatz 
for the boson and fermion density distributions. We will 
assume that the radial profile of the fermion and boson 
gases is a gaussian. This is an exact result in the quasi- ID 
limit (formally expressed as hio^ ^ /^b and Tiloy- S> /^f, 
where /^b and /ip are the chemical potentials of the boson 
gas and fermion gas, respectively |32j ) for which the ra- 
dial profile coincides with the gaussian ground harmonic 
oscillator state. In the following, we do not impose a 
fixed size to our gaussian profile but simply require the 
transverse profile to approximate a gaussian shape. Thus 
our approach can apply to systems that lie outside of this 
dimensionality condition. However, as is shown later, the 
radial sizes of the solitons remain close to the radial har- 
monic oscillator length Zho throughout. 

The most obvious choice for a suitable axial profile is, 
by analogy to ID bright bosonic soliton result, a sech- 
profile [H]. However, with this choice we are unable 
to obtain analytic solutions for the variational energies. 
Karpiuk et al. |18j pursued this choice numerically. In- 
stead, to obtain an analytic form for the variational en- 
ergies, we employ a gaussian axial profile. From stud- 
ies of bright BEC solitons it has been shown, firstly. 



that sech and gaussian axial profiles give very similar re- 
sults, and secondly, that both forms of ansatz give very 
good agreement with more precise theoretical treatments, 
e.g. numerical solutions of the Gross-Pitaevskii equa- 
tion [TU] [TT] . We will thus consider the boson {fermion} 
density nB{F} to have a cylindrically-symmetric gaussian 
profile, 

-B{F}(r) = ^§g^exp(-ycxp(-^), (1) 

where Lr and Lz are the radial and axial sizes, respec- 
tively, and A^B{F} is the numbers of bosons {fermions}. 
We consider that Lr and Lz are common for both the 
bosons and the fermions. 

Note that the validity of the mean-field Gross- 
Pitaevskii model for the boson gas requires that A'b ^ 1. 
Furthermore, our description of the Fermi gas component 
is based on the Thomas-Fermi approximation which is 
valid for TVp' 1 [21]. Hence our approach is strictly 
valid only in the regime of large N and we will only con- 
sider parameters that satisfy this. 

B. Energetics of the system 

We will consider the total energy density of the Bose- 
Fermi state e[nB,nF] to be the sum of the boson con- 
tribution eB["-B], fermion contribution £p[np] and boson- 
fermion term ebf [^b , "-f] |23) . We will proceed by mod- 
elling each energy contribution of the gaussian wavepack- 
ets in turn. Note that the energy is the volume integral 
of the corresponding energy density E = J e[n]dV. For 
a different choice of ansatz and in the absence of p-wave 
interactions, Karpiuk et al. |18j followed a similar varia- 
tional approach to explore Bose-Fermi soliton solutions. 

L Bosonic energy contribution 

The energy density of the zero-temperature boson gas, 
interacting via s-wave interactions of scattering length 
Ob, is provided by the Gross-Pitaevksii model |2F, 

£bM - ^|V0iip + l/B(r)nB + ^^^4. (2) 

The terms of the right-hand side represent, respectively, 
the kinetic, potential and interaction energies of the bo- 
son gas. For convenience we will express length in terms 
of the boson harmonic oscillator length l^o — \/h/mBUJB 
and adopt dimensionless variational lengthscales lz = 
Lz/lho and lr — Lr/lho- Furthermore, we rescale en- 
ergy by the bosonic harmonic oscillator energy fkoB via 
E — >• E/{hujB)- Substituting nB(r) into Eq. ^ and inte- 
grating over space, one arrives at the well-known expres- 
sion for the total boson energy [21], 

+ + (3) 

A^B 2 \Pr ^IV 2 ^ V2nko Ur 
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From this the presence of the collapse instability for at- 
tractively interacting BECs (as < 0) is clear: the inter- 
action term (negative) dominates over all other (positive) 
energy contributions as the size of the wavepacket tends 
towards zero. While meta-stable states can be formed for 
certain parameters, the lowest energy state described by 
this equation is always a collapsed state of zero width. 



2,. Fermionic energy contribution 



where aeps is the s-wave boson-fermion scattering co- 
efficient, flBFp is the p-wave scattering length and jjL = 
niBmp /{rriB +mp ) . For the density profiles ([T]) the boson- 
fermion interaction energy is. 
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Determination of the mean-field fermionic energy den- 
sity in general requires solving N-p coupled Hartree-Fock 
equations. In this manner, Karpiuk et al. |18| success- 
fully modelled BF soliton solutions, but the approach 
is only tractable for of the order of 10 fcrmions. In- 
stead, we will adopt an analytic form for the fermion 
energy density derived by Roth and Feldmeier ^23]. Em- 
ploying the Thomas-Fermi approximation, this approach 
follows from deriving the energy density of a homoge- 
neous fermionic system and replacing the fixed density 
with nY{v) (thus neglecting the contribution to the en- 
ergy from density gradients) . The Thomas-Fermi approx- 
imation is valid for the equlibrium state of a fermi gas 
when the fermion wavelength is much smaller than the 
system size. For a trapped fermi gas, it can be shown 

1 /3 

that this is satisfied when N-p ^ 1 [5T]. In agree- 
ment with this, the Thomas- Fermi approximation has 
been shown to give an excellent mean-field description of 
BF mixtures (in close agreement with Hartree-Fock cal- 
culations) for A^F ^ 1000 [23l |33]. Via this approach, 
the energy density of the fermions (potential, kinetic and 
p-wave interaction terms, respectively) is 23 , 
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where uf is the p-wave contact interaction between 
fermions [23]. Inserting the gaussian profile riF(r) and 
integrating gives the total fermionic energy, 
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where a = (9/50)62/3(3/5)^/2^^/3 w 0.6743 and /3 
(3/160^4)(3^/8)i/2(g^2)5/3 ^ 0.1880. 



3. Boae-Fermi energy contribution 

Again under the Thomas-Fermi approximation, the in- 
teraction energy density between the bosons and fermions 
eBFlnBi'rap] is [35], 



eBF[nB, ?if] 
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Energy landscapes and obtaining the variational 
solutions 



The total energy of our gaussian Bose-Fermi wavepack- 
ets is given hy E = Eb + Ef + Ebf- Upon fixing the 
experimental parameters (atom masses, atom numbers, 
scattering lengths and trap frequencies), the energy be- 
comes confined to being a function of only the wavepacket 
size parameters Ir and Iz- This function E{lr,lz) can be 
visualised as an energy landscape in which the presence 
of an energy minimum represents a variational solution. 
In practice, we numerically define this energy landscape 
and perform a simple computational search for such en- 
ergy minima. Note that an unphysical solution can oc- 
cur at the origin representing the effect of collapse. It 
is unphysical in the sense that a real Bose-Fermi system 
cannot shrink to zero size; in reality, a collapse will even- 
tually become halted by the surge of three-body losses 
as the gas densities rise. At most, only one physical so- 
lution is ever present in these energy landscapes. We 
denote the coordinates of such a variational solution by 
the coordinates and 

Where we map out regions of soliton solutions within 
a particular parameter space, e.g., obfs — Nf space, this 
is done by randomly sampling combinations of these pa- 
rameters. We typically restrict our numerical search to 
landscapes of extent [0, 2] Ir x [0, 10] Iz- 

We have verified that our methodology gives very close 
agreement with the variational stability diagram pre- 
sented in Fig. 7 of [TB]. Our results deviate by less than 
10%, a deviation which is to be expected given the dif- 
ferent ansatz employed. 



D. Analytical limit: no p-wave interactions and 

Nb > iVp 

We can gain a simplified analytic form for the total 
variational energy if we neglect p-wave interactions (of = 
OBFp = 0) and assume A'^b ^ -^f- The latter condition 
renders the fermion-fermion energy terms negligible and 
makes significant only terms involving A^b- The total 
variational energy then reduces to. 
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This form wiU enable us to gain physical intuition of 
the system and provide simple criteria for the existence 
of stable Bose-Fermi solitons. Importantly, it has the 
same form as the gaussian variational energy of a purely 
bosonic gas [lOj, but with an effective s-wave scattering 
length given by [24], 



Oeff = Ob + flBFs 



iVp 



(9) 



A rudimentary requirement for the ability to self-trap is 
that the net interactions are attractive {aos < 0). This 
places a lower bound on the ratio Np/N-g for which Bose- 
Fermi solitons can be self-supported. 
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Furthermore, bright bosonic solitons are established to 
collapse when the scattering length is less than the criti- 
cal value Ob = fcc'ho/-^B- For a 3D gaussian wavepacket 
the dimensionless coefficient has been shown to be kc = 
—0.778 [inj (for other shapes of wavepacket the value dif- 
fers but remains of the order of unity) . This leads to an 
upper bound for the ratio Np/Ns in order to prevent 
collapse of the system. 
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From Eqs. ( 10|11 ) it is evident that the soliton solu- 
tions exist within a "window" of fermion atom number 
Np whose width is. 



ANp 



ho 



a-BFs 



mp 



mp + mB 
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Equations (10), (11) and (12) provide us with an es- 



timate for the locality and range over which BF solitons 
solutions may exist. Equation (12) indicates that the 



width of the soliton bands can be extended by employing 
weaker radial trapping and weaker Bose-Fermi scattering 
length. However, we will find that the soliton bands are 
even more restricted if the condition Nb/Np ^ 1 is re- 
moved. Indeed, the above equations predict the existence 
of soliton solutions at the native Bose-Fermi scattering 
length ubfs — — 215ao. However, as we shall see, the full 
variational results predict that the soliton bands become 
vanishing narrow at this scattering length. 



III. RESULTS 

While our analytical results presented so far are appli- 
cable to any Bose-Fermi species, we will focus our ensuing 
results on a *^Rb and ""^K mixture due to its naturally 
strong and attractive Bose-Fermi coupling and its promi- 
nent experimental occurence to date [T31 [IS]. We will 
assume that the atoms are spin-polarized and confined 



(in the radial direction) by a magnetic trap such that the 
trap frequencies are related via wf/wb = (wb/tif)^^^. 
Throughout our results we fix the boson-boson s-wave 
scattering length to be the experimentally measured 



value of Ob — 99ao 



where gq — 5.3 x 10 m is 



the Bohr radius. We will consider the radial trap fre- 
quency, boson and fermion atom numbers, and remaining 
scattering lengths to be variables. Note that scattering 
lengths can typically be experimentally varied over many 
orders of magnitude by using magnetic and optical fields 
to access inter-atomic scattering resonances. In the cases 
where we fix the boson-fermion s-wave scattering length, 
we take it to be its native value of gbfs = — 215ao, as 
measured by Ferlaino et al. 31 . 

First we will explore the soliton solutions in the ab- 
sence of p-wave interactions. Following this we will con- 
sider how p-wsye interactions between fermions modify 
the soliton solutions. We will not explicitly present re- 
sults for finite boson-fermion p-wave interactions but will 
comment on how they affect the system in Section IV. 



A. Absence of p-wave interactions 

1. Soliton hands in N-p — sbps space 

Figure [ija) presents the soliton solutions in the pa- 
rameter space of aBFs — ^f for 10^ bosons and various 
radial trap frequencies [a;B/27r — 10 (green solid region), 
100 (blue dotted region) and 1000 Hz (black dashed re- 
gion)]. For each trap frequency, the soliton solutions ex- 
ist in a narrow band in the obfs < half-plane. (Note 
that on the scale of the figure these bands appear as 
lines). Above each band, the system is unstable to dis- 
persion and below the band it is unstable to collapse. 
Each band scales approximately as —1/Np, as indicated 
by rearranging Eqs. ( 10 ) and ( 11 ). The bands are weakly 
dependent on trap frequency. Indeed their differences 



are only apparent on a much more magnified scale (in- 
set of Figjlja)). As ujb is increased, the bands shift up- 
wards in gbps and become slightly narrower. The latter 
change is in qualitative agreement with equations ( |10p2 ) 
which predict the band width to scale in proportion to 
Zho — \/fi/mujB- Although not visible in Figure [l](a), the 
bands become progressively narrower as Np increases. 
Indeed, beyond some critical fermion number Np"^ we 
cannot detect further solutions. This marks an important 
difference between the approximated analytic predictions 
of Eqs. (10)-(12| and the full variational solutions. For 



example, for trap frequencies ujB/2n — 10, 100 and 1000, 
A^^"* « 760, 660 and 620, respectively. 

Within the soliton bands, the soliton radial size re- 
mains close to ^ho throughout. The axial size l'^ is infinite 
at its dispersive boundary and reduces as gbps is made 
more attractive, with the soliton becoming almost spher- 
ical at the point of collapse. This is qualitatively similar 
to the case for BEC bright solitons [TU] . 

Figure [2] shows how the soliton bands change for differ- 
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FIG. 1: (a) The bands of soliton solutions in aeps — A^f space 
for a *'^Rb-*''K mixture with zero p-wave interaction. We 
present a fixed number of bosons A^B = 10^ and various trap 
frequencies tJB/27r =10 (green sohd), 100 (blue dotted), 1000 
(black dashed) Hz. The inset provides a close up of the soliton 
bands, (b) Energy landscapes {ijj-b/2-k = lOHz, A^f = 500 
and A''b = W) showing four distinct regimes: (i) well above 
the soliton band (obfs = — 6200ao), (ii) just above the band 
(ciBFs = — 6260ao), (iii) within the band (obfs = — 6296ao) 
and (iv) below the band, (sbfs = — 6350ao). These plots 
corresponds to the crosses in (a). In (iii) the soliton solution 
is highlighted by the white cross. 



bounds for soliton-supporting conditions to be met, see, 
e.g., Eqs. (10l-(12)) is decreased the soliton solutions ex- 
tend to smallerJaeFs I ■ Indeed, if one were to extrap- 
olate our predictions (beyond its strict regime of valid- 
ity) to lower atom number, one could imagine that the 
soliton bands would reach obps — — 215ao. Indeed, us- 
ing a model valid for low atom numbers, Karpiuk et al. 
[IH] predict the existence of BF solitons at the native 
obfs — — 215ao for very low atom numbers. 



2. Energy landscape regimes 

To gain physical insight into the system, in Fig.[T]^b)(i)- 
(iv) we present energy landscapes of four distinct regimes 
in this parameter space. While we present the land- 
scapes for a specific set of parameters (wb/Stt ~ 10 Hz, 
N-p = 500 and N-q = 10^) the qualitative behaviour is 
generic. The location of each case is indicated in the inset 
of Fig.jIJa) by crosses. These four regimes [Fig.[T]^i)-(iv)] 
are as follows: 

• (i) Sufficiently above the soliton band, the net con- 
tact interactions are repulsive (positive) and the en- 
ergy landscape consists of a downward 'chute' aligned 
along the axis. Any wavepacket subjected to this 
system will disperse axially. 

• (ii) Just above the soliton band the net interactions be- 
come attractive (negative) and compete with the pos- 
itive energy terms. The chute remains but the global 
energy minimum is now at the origin, where the energy 
diverges to -co. The attraction is insufficiently strong 
to support a soliton. 

• (iii) Within the band, the play-off between the inter- 
species attraction and repulsive zero-point kinetic en- 
ergy leads to a local energy minimum at [/°, l^] (high- 
lighted by the white cross), corresponding to the self- 
trapped soliton solution. 

• (iv) Below the soliton band the attractive interactions 
dominate the kinetic energy such that the only energy 
minimum is the global minimum at the origin, repre- 
senting the collapse instability of the system. 



ent number of bosons (in both linear (a) and logarithmic 
plots (b)). For increasing N-q the bands shift to more neg- 
ative flBFs and larger N-p, and the bands becomes wider 
as A^B is reduced. 

Our numerical results predict that, for the numbers 
of bosons and fermions permitted by our model (3> 1), 
soliton solutions do not occur at the native Bose-Fermi 
scattering length asPs — — 215ao. The soliton solutions 
exist only for scattering lengths obfs ^ — 215ao. How- 
ever, our results show that as A'b, and by association 
A^F (since the ratio N-q/Np must remain within narrow 



3. Comparison to the simplified limit of Eqs. 10) and ,11) 



In Fig.[2jb) we present a comparison between the varia- 
tional predictions for the soliton bands and the simplified 



analytic estimates provided by Eqs. ( 10 ) (dashed line) 



and (11) (dotted line). For ease of observing the differ- 
ences between the two methods, the data is presented on 
a log- log plot. For low number of bosons (A''^ = 500) the 
variational solutions deviate from the predictions. How- 
ever, for larger boson number the agreement is excellent. 
This is to be expected since these predictions assume 



6 




(a) 200 



1000 



FIG. 2: (a) The bands of soliton solutions in aeps — Np 
space for no p-wave contributions. We consider fixed trapping 
cjb/Ztt = lOOHz and various boson numbers of A^b = 500 
(green), 5,000 (black), 50,000 (red) and 100,000 (blue), (b) 
Log-log plot of the sol iton bands in (a), with the analytic 
predictions of Eqs. (10 1 (dashed line) and (111 (dotted line). 



The grey arrows indicate the direction of increasing A^b. 



Nb/N-p ^ 1. Indeed for Nb = 5, 000 and above, the sim- 
plified forms give very good predictions for the regimes 
of soliton solutions (indistinguishable from the full vari- 
ational results on the scale of this figure). 

An important difference, however, is that according 
to the analytical result, the width of the soliton band 
decreases as 1/Np, becoming vanishingly small only as 
A^F ^ oo. In contrast, the full variational solutions dis- 
appear beyond a finite Np. Furthermore, numerically, 
the bands increase in width as N-q is increased, whereas 
the analytical prediction for the band width is indepen- 
dent of Nn- 
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FIG. 3: (a) Soliton bands in ap — Nf space for boson- 
fermion interaction asps ~ — 215ao in presence of fermion- 
fermion p-wave interactions. The number of bosons is fixed 
to A^B = 100, 000 and we present trap frequencies of i^b/Stt = 
10 (horizontal thatch), 100 (vertical thatch) and 1,000 Hz 
(grey region), (b) Energy landscapes (for ujb/2tt = 100 and 
Nf = 20,000) of (i) below the band (ap = 50ao) (ii) within 
the band (ap = 125ao) and (iii) above the band (ap — 180ao). 
The locations of these cases are indicated in plot (a). 



B. The role of p-wave fermion interactions 

1. Soliton bands in ap — Np space 

We now consider the presence of fermion-fermion (p- 
wave) interactions. We fix the boson- fermion s-wave in- 
teraction to its natural value a^Fs = — 215ao and explore 
the parameter space of ap — N-p . The results are shown in 
Fig.wlfor fixed boson number Nb ~ 10^ and various trap 
frequencies [wb/Ztt = 10 (horizontal thatch), 100 (verti- 
cal thatch) and 1000 Hz (grey region)]. Recall that in 
the absence of p-wave interactions, no solitons were ob- 
tained for flBFs = — 215ao. In contrast, in the presence 
of repulsive p-wave fermion-fermion interactions, we now 
see extensive regions of soliton solutions. The regions be- 
come larger for increased trap frequency. This change in 
size occurs due to a shift in the lower boundary of the 
regions; the upper boundary is insensitive to ujb, as can 
be seen in Fig. [Sj^a). 

This enhanced stability when p-wave fermion-fermion 
interactions are included arises from the fact that the 
fermion-fermion interaction term in Eq. ([5| scales as 
Thus, if ap > 0, this term ensures that the en- 
ergy diverges to positive values as / — > and completely 



removes the presence of a collapse instability. This ener- 
getic behaviour is demonstrated by the landscapes shown 
in Fig.^b). 

• (i) Below the relevant soliton band (e.g. the band 
shaded with vertical hatch in Fig. [sj^a)) a significant 
collapse region is present and dominates the energy 
landscape. However, a positive region is just visible 
close to the origin in Fig. [3][b)(i). 

• (ii) Solutions become supported when the fermion- 
fermion interaction becomes larger and the play-off be- 
tween all of the energy contributions generates a local 
minimum in the energy landscape (white cross). 

• (iii) Above the soliton band, the repulsive fermion- 
fermion interaction becomes so large that it makes the 
system fully dispersive. 

Case (ii) is an intriguing prediction. It suggests that 
the presence of repulsive p-wave fermion interactions 
leads to solitons which are the global energy minimum 
of the system. This indicates that such solitons would 
be far more robust and stable than their bosonic coun- 
terparts, which are well-known to exist as meta-stable 
states prone to an irremovable collapse instability. 
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FIG. 4: Soliton bands in op — Np space for fixed trap fre- 
quency oje/Stt = 100 Hz and boson numbers of A''b — 10, 000 
(solid), 50,000 (dashed) and 100,000 (dotted). 



In Fig.|4] we show how the bands change with the num- 
ber of bosons. As the number of bosons decreases, the 
bands shift to lower Np and become narrower. Indeed, 
the bands scale approximately as 1/Nb, i.e. if we plot 
Nf/Nb on the x-axis, the bands approximately overlap 
with each other. 



C. Soliton bands in ap — hbps space 

The simultaneous manipulation of more than one scat- 
tering length has not been experimentally demonstrated. 
However, in principle, this could be possible through a 
combination of magnetic, optical and confinement reso- 
nances. With this is mind and by way of exploring the 
soliton solutions further, we now turn to examine the re- 
gions of soliton solutions in ap — obfs space [Fig. [sjja)], 
for fixed A'b and Np. We will see that parameter space 
is particularly interesting because it offers the possibility 
of forming bright soliton solutions as clear global minima. 
In this case, we observe complex-shaped regions of soli- 
ton solutions. This includes a large cusp shaped region in 
the ap > half plane. Although not shown, this region 
extends to indefinitely negative obfs- Furthermore, the 
soliton region features narrow 'fingers' which extend far 
into the ap < half-plane. 

We discriminate six distinct regions in this parameter 
space which we interpret by reference to their typical 
energy landscapes presented in Fig. [5|b): 

• (i) Above this narrow band the landscape is dominated 
by dispersion with a localised collapse region, but no 
local minimum exists. 

• (ii) Within the soliton band there exists a shallow en- 
ergy minimum (case(ii)) adjacent to the collapse and 
dispersive regions. 
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FIG. 5; Soliton bands (as shown by their boundary lines) in 
the parameter space of ap and obps. (a) We fix A''b = 100, 000 
and A'p — 500, and use ^jb/Ztt — 10 (red dashed), 100 (green 
solid) and 1, 000 (blue dot-dashed) Hz. (b) Energy landscapes 
for the wb/Stt = 10 case at the points (i) ap = — 2000ao, 
"bps = — 6280oo, (ii) ap = — 2000ao, asps ~ — 6295ao, (iii) 
ap = — 2000ao, aBPs = — 6320ao, (iv) ap = 2000ao, asps = 
— 6280oo, (v) ap = SOOOao, aaps = — 6320ao, and (iv) ap = 
4000ao, aBPs = — 6280ao. These points are denoted in (a). 



• (iii) Below the soliton band the whole landscape is un- 
stable to collapse. 

• (iv, v) In the regions containing points (iv) and (v) 
there is no collapse region at the origin and there exists 
a well-localised and deep energy minimum denoting a 
soliton solution. 

• (vi) This region is dispersive due to the dominance of 
repulsive interactions. 

Regions (i),(ii),(iii) and (vi) possess energy landscapes 
which are analogous to those seen in Fig. [ijb). How- 
ever, the most intriguing regions are (iv) and (v). These 
solutions are the most common type that exist in this 
parameter space. Like the observation in Fig. [3]jb)(ii), 
the soliton now becomes the global energy minimum of 
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FIG. 6: Soliton bands in the parameter space ap — obps 
(shown by their boundary hnes). We fix ojb/Stt = lOOHz and 
iVp = 500 and consider Nb = 1,000 (dashed fine), 10,000 
(sohd line) and 100,000 (dotted line). 



the system. However, these landscapes are strikingly well 
localised and deep. The depth of this minimum is typi- 
cally of the order of lOOhu^ . For comparison, for a bright 
bosonic soliton the depth of the energy minima is of the 
order of OAhuj^- Within the context of our scaling solu- 
tions (fixed gaussian shape), this depth and narrowness 
of the energy minima indicates extreme stability of the 
solutions to shape modification, including collapse and 
dispersion. 

For completenesss, Fig. [6] demonstrates how these ex- 
tensive soliton regions in ap — a^Ps space become mod- 
ified for different Nb- The presence of the "fingers" is 
sensitive to but the main region of solutions persists, 
albeit shifting to more negative a^Fs with increasing A^b • 

While here we have limited our study to the p-wave 
interactions of only the fermions, the same qualitative 
soliton regimes and solutions are obtained if the p-wave 
boson-fermion interaction is included instead (or, indeed, 
if both are included) . This is because the p-wave fermion- 
boson interaction has the same functional form (scaling 

as r^). 



IV. DISCUSSION AND CONCLUSION 

Our results show that in the absence of p-wave con- 
tributions the range of soliton solutions is rather lim- 
ited, being confined to narrow bands in obf ~ -^f space. 
Indeed, the solutions behave similarly to bright bosonic 
solitons but with an effective scattering length. While re- 
ducing the radial confinement widens the soliton bands, 
this also makes the system more prone to collective ex- 
citations which may disrupt the soliton. For a ^^Rb- 
mixture and in the validity regime of our approach 
(A^B,-^F ^ 1) we cannot locate soliton solutions at the 
native boson-fermion scattering length of — 215ao. Values 



of around — lOOOao are required to reach a soliton band at 
these atom numbers. The work of Karpiuk et al. demon- 
strated that atom numbers of order unity are required 
to support Bose-Fermi solitons at the native ^''Rb-'*°K 
scattering length. While the scattering length can be 
engineered to large values through scattering resonance 
tuning, the ensuing soliton bands are so narrow that ex- 
perimentally locating them is likely to be problematic. 
Furthermore, there remains a collapse instability in the 
system and the ratio of bosons to fermions is constrained 
to small values. In contrast, the presence of repulsive 
p-wave fermion-fcrmion interactions has a dramatic sta- 
bilizing effect on the system. This can lead to a removal 
of the collapse instability such that the soliton solutions 
become global energy minima. We find extensive soliton 
regimes, in which the soliton minima are extremely deep, 
suggesting that they may form soliton structures that 
are considerably more robust than in the absence of p- 
wave interactions. The p-wave interaction also provides a 
strong tuning parameter, enabling the boson-fermion ra- 
tio to be dramatically varied. Importantly, for a potential 
experimental realization of BF solitons, we find extensive 
soliton solutions at the native boson-fermion interaction 
and with only moderate fermionic interactions. 

The remarkable capacity of repulsive p-wave interac- 
tions to remove the collapse instability stems from the 
scaling behaviour of its energy contribution. Denoting a 
generalized lengthscale of the gaussian wavepacket as 
the total variational energy of the Bose-Fermi system is 
of the form, 



E 



1 



1 



1 



(13) 



The first two terms, the kinetic and potential energy 
terms, are always positive. The last two terms, the s- 
wave and p-wave interaction energies, respectively, may 
be positive or negative. In the absence of p-wave inter- 
actions, a negative s-wave term will cause the energy to 
diverge to — oo as £ — >■ 0, signifying the presence of the 
collapse instability. For non-zero p-wave interactions, the 
p-wave term dictates the fate of the system as £ — > and 
importantly, for positive p-wave interactions the collapse 
divergence is completely removed. 

It is important to note that this scaling behaviour orig- 
inates from the Thomas-Fermi approximation and so is 
not limited to gaussian wavepackets. Consider homoge- 
neous Bose and Fermi gases in a large hard-wall box of 
volume i^. It is trivial to see from Eqs. S and (|4| that 
the energy scales as above, minus the ^^erm. Thus it 
is clear that the capacity of repulsive p-wave interactions 
to stabilise against collapse is likely to extend to trapped 
Bose-Fermi mixtures in general. 

While we have presented results for p-wave interac- 
tions in only the fermion-fermion case, we find quali- 
tatively similar soliton regions, landscapes and conclu- 
sions when including boson-fermion p-wave interactions 
instead. This is because the same energy scaling dis- 
cussed above applies. 
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In conclusion, according to a variational model, the 
presence of repulsive p-wave interactions in Bosc-Fcrmi 
mixtures removes the usual collapse instability and leads 
to stable, robust bright soliton solutions that are global 
energy minima of the system. We have discussed specifi- 
cally the boson-fermion pairing of ^''Rb and ^°K, but the 
stabilizing effect of repulsive p-wave interactions should 



apply in general to other ultracold Bosc-Fermi mixtures. 
Given that the collapse instability has proved a major 
hindrance to the controlled generation, manipulation and 
interaction of matter-wave solitons to date, these more 
stable p-wme entities may provide a more versatile route 
to explore and exploit the special characteristics of soli- 
tons. 
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